home *** CD-ROM | disk | FTP | other *** search
/ Aminet 28 / Aminet 28 (1998)(GTI - Schatztruhe)[!][Dec 1998].iso / Aminet / dev / lang / fpc09905c.lha / fpc / inc / math.inc < prev    next >
Text File  |  1998-09-21  |  34KB  |  943 lines
  1. {
  2.     $Id: math.inc,v 1.1.1.1 1998/03/25 11:18:44 root Exp $
  3.     This file is part of the Free Pascal run time library.
  4.     Copyright (c) 1993,97 by several people
  5.     member of the Free Pascal development team.
  6.  
  7.     See the file COPYING.FPC, included in this distribution,
  8.     for details about the copyright.
  9.  
  10.     This program is distributed in the hope that it will be useful,
  11.     but WITHOUT ANY WARRANTY; without even the implied warranty of
  12.     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
  13.  
  14.  **********************************************************************}
  15. {*************************************************************************}
  16. {  math.inc                                                               }
  17. {*************************************************************************}
  18. {       Copyright Abandoned, 1987, Fred Fish                              }
  19. {                                                                         }
  20. {       This previously copyrighted work has been placed into the         }
  21. {       public domain by the author (Fred Fish) and may be freely used    }
  22. {       for any purpose, private or commercial.  I would appreciate       }
  23. {       it, as a courtesy, if this notice is left in all copies and       }
  24. {       derivative works.  Thank you, and enjoy...                        }
  25. {                                                                         }
  26. {       The author makes no warranty of any kind with respect to this     }
  27. {       product and explicitly disclaims any implied warranties of        }
  28. {       merchantability or fitness for any particular purpose.            }
  29. {-------------------------------------------------------------------------}
  30. {       Copyright (c) 1992 Odent Jean Philippe                            }
  31. {                                                                         }
  32. {       The source can be modified as long as my name appears and some    }
  33. {       notes explaining the modifications done are included in the file. }
  34. {-------------------------------------------------------------------------}
  35. {       Copyright (c) 1997 Carl Eric Codere                               }
  36. {                                                                         }
  37. {*************************************************************************}
  38. {  This is the Motorola 680x0 specific port of the math include.          }
  39. {*************************************************************************}
  40. {                                                                         }
  41. {  o all reals are mapped to the single type under the motorola version   }
  42. {                                                                         }
  43. {   What is left to do:                                                   }
  44. {    o add support for sqrt with fixed.                                   }
  45.  
  46. type
  47.     TabCoef = array[0..6] of Real;
  48.  
  49.  
  50. const
  51.       PIO2   =  1.57079632679489661923;       {  pi/2        }
  52.       PIO4   =  7.85398163397448309616E-1;    {  pi/4        }
  53.       SQRT2  =  1.41421356237309504880;       {  sqrt(2)     }
  54.       SQRTH  =  7.07106781186547524401E-1;    {  sqrt(2)/2   }
  55.       LOG2E  =  1.4426950408889634073599;     {  1/log(2)    }
  56.       SQ2OPI =  7.9788456080286535587989E-1;  {  sqrt( 2/pi )}
  57.       LOGE2  =  6.93147180559945309417E-1;    {  log(2)      }
  58.       LOGSQ2 =  3.46573590279972654709E-1;    {  log(2)/2    }
  59.       THPIO4 =  2.35619449019234492885;       {  3*pi/4      }
  60.       TWOOPI =  6.36619772367581343075535E-1; {  2/pi        }
  61.       lossth =  1.073741824e9;
  62.       MAXLOG =  8.8029691931113054295988E1;    { log(2**127)  }
  63.       MINLOG = -8.872283911167299960540E1;     { log(2**-128) }
  64.  
  65.       DP1 =   7.85398125648498535156E-1;
  66.       DP2 =   3.77489470793079817668E-8;
  67.       DP3 =   2.69515142907905952645E-15;
  68.  
  69. const sincof : TabCoef = (
  70.                 1.58962301576546568060E-10,
  71.                -2.50507477628578072866E-8,
  72.                 2.75573136213857245213E-6,
  73.                -1.98412698295895385996E-4,
  74.                 8.33333333332211858878E-3,
  75.                -1.66666666666666307295E-1, 0);
  76.       coscof : TabCoef = (
  77.                -1.13585365213876817300E-11,
  78.                 2.08757008419747316778E-9,
  79.                -2.75573141792967388112E-7,
  80.                 2.48015872888517045348E-5,
  81.                -1.38888888888730564116E-3,
  82.                 4.16666666666665929218E-2, 0);
  83.  
  84.  
  85.  
  86.  
  87.     function int(d : real) : real;
  88.       begin
  89.         { this will be correct since real = single in the case of }
  90.         { the motorola version of the compiler...                 }
  91.         int:=real(trunc(d));
  92.       end;
  93.  
  94.     function trunc(d : real) : longint;
  95.     var
  96.      l: longint;
  97.     Begin
  98.      asm
  99.         move.l   d,d0           { get number                        }
  100.         move.l   d2,-(sp)       { save register                     }
  101.         move.l   d0,d1
  102.         swap     d1             { extract exp                       }
  103.         move.w   d1,d2          { extract sign                      }
  104.         bclr     #15,d1         { kill sign bit                     }
  105.         lsr.w    #7,d1
  106.  
  107.         and.l    #$7fffff,d0    { remove exponent from mantissa     }
  108.         bset     #23,d0         { restore implied leading "1"       }
  109.  
  110.         cmp.w    #BIAS4,d1      { check exponent                    }
  111.         blt      @zero           { strictly factional, no integer part ?   }
  112.         cmp.w    #BIAS4+32,d1   { is it too big to fit in a 32-bit integer ? }
  113.         bgt      @toobig
  114.  
  115.         sub.w    #BIAS4+24,d1   { adjust exponent                   }
  116.         bgt      @trunclab2     { shift up                          }
  117.         beq      @trunclab7     { no shift (never too big)          }
  118.  
  119.         neg.w    d1
  120.         lsr.l    d1,d0          { shift down to align radix point;  }
  121.               { extra bits fall off the end (no rounding) }
  122.         bra      @trunclab7      { never too big                     }
  123.     @trunclab2:
  124.         lsl.l   d1,d0           { shift up to align radix point     }
  125.     @trunclab3:
  126.         cmp.l   #$80000000,d0   { -2147483648 is a nasty evil special case }
  127.         bne      @trunclab6
  128.         tst.w    d2             { this had better be -2^31 and not 2^31    }
  129.         bpl      @toobig
  130.         bra      @trunclab8
  131.     @trunclab6:
  132.         tst.l   d0              { sign bit set ? (i.e. too big)     }
  133.         bmi     @toobig
  134.     @trunclab7:
  135.         tst.w   d2              { is it negative ?                  }
  136.         bpl     @trunclab8
  137.         neg.l   d0              { negate                            }
  138.         bra     @trunclab8
  139.     @zero:
  140.         clr.l   d0              { make the whole thing zero         }
  141.         bra     @trunclab8
  142.     @toobig:
  143.         moveq   #-1,d0          { ugh. Should cause a trap here.    }
  144.         bclr    #31,d0          { make it #0x7fffffff               }
  145.     @trunclab8:
  146.         move.l  (sp)+,d2
  147.         move.l  d0,l
  148.      end;
  149.      if l = $7fffffff then
  150.       RunError(207)
  151.      else
  152.        trunc := l
  153.     end;
  154.  
  155.  
  156.  
  157.  
  158.  
  159.     function abs(d : Real) : Real;
  160.     begin
  161.        if( d < 0.0 ) then
  162.          abs := -d
  163.       else
  164.          abs := d ;
  165.     end;
  166.  
  167.  
  168.     function frexp(x:Real; var e:Integer ):Real;
  169.     {*  frexp() extracts the exponent from x.  It returns an integer     *}
  170.     {*  power of two to expnt and the significand between 0.5 and 1      *}
  171.     {*  to y.  Thus  x = y * 2**expn.                                    *}
  172.     begin
  173.       e :=0;
  174.       if (abs(x)<0.5) then
  175.        While (abs(x)<0.5) do
  176.        begin
  177.          x := x*2;
  178.          Dec(e);
  179.        end
  180.       else
  181.        While (abs(x)>1) do
  182.        begin
  183.          x := x/2;
  184.          Inc(e);
  185.        end;
  186.       frexp := x;
  187.     end;
  188.  
  189.  
  190.     function ldexp( x: Real; N: Integer):Real;
  191.     {* ldexp() multiplies x by 2**n.                                    *}
  192.     var r : Real;
  193.     begin
  194.       R := 1;
  195.       if N>0 then
  196.          while N>0 do
  197.          begin
  198.             R:=R*2;
  199.             Dec(N);
  200.          end
  201.       else
  202.         while N<0 do
  203.         begin
  204.            R:=R/2;
  205.            Inc(N);
  206.         end;
  207.       ldexp := x * R;
  208.     end;
  209.  
  210.  
  211.     function polevl(var x:Real; var Coef:TabCoef; N:Integer):Real;
  212.     {***************************